# 基于《基于OpenSees的某多层钢框架结构的抗震分析_李文海.pdf》
# 4层3×2跨钢框架结构 - 国际单位制版本（m-N-kg）

wipe
model basic -ndm 3 -ndf 6

# 材料定义（严格按文档参数，国际单位制：m-N）
# 混凝土材料 Concrete02（单位：N/m²）
# fc'=26.8MPa=26.8×10⁶N/m², epsc0=0.002, fcu'=10MPa=10×10⁶N/m², epsu=0.006, lambda=0.1, ft=2.2MPa=2.2×10⁶N/m², Ets=44000MPa=44×10⁹N/m²
uniaxialMaterial Concrete02 1 -26.8e6 -0.002 -10.0e6 -0.006 0.1 2.2e6 44.0e9

# 钢材材料 Steel01 (Q345钢)（单位：N/m²）
# Fy=345MPa=345×10⁶N/m², E0=206000MPa=206×10⁹N/m², b=0.01
uniaxialMaterial Steel01 2 345.0e6 206.0e9 0.01

# 几何参数（国际单位制：米）
set L1 4.0;     # X向跨度 4m
set L2 3.5;     # Y向跨度 3.5m  
set H  3.0;      # 层高 3m

# 节点创建 - 4层3×2跨框架
# 底层节点 (Z=0)
for {set i 0} {$i <= 3} {incr i} {
    for {set j 0} {$j <= 2} {incr j} {
        set nodeTag [expr $i*3 + $j + 1]
        node $nodeTag [expr $i*$L1] [expr $j*$L2] 0.0
    }
}

# 二层节点 (Z=3.0)
for {set i 0} {$i <= 3} {incr i} {
    for {set j 0} {$j <= 2} {incr j} {
        set nodeTag [expr $i*3 + $j + 13]
        node $nodeTag [expr $i*$L1] [expr $j*$L2] [expr 1*$H]
    }
}

# 三层节点 (Z=6.0)
for {set i 0} {$i <= 3} {incr i} {
    for {set j 0} {$j <= 2} {incr j} {
        set nodeTag [expr $i*3 + $j + 25]
        node $nodeTag [expr $i*$L1] [expr $j*$L2] [expr 2*$H]
    }
}

# 四层节点 (Z=9.0)
for {set i 0} {$i <= 3} {incr i} {
    for {set j 0} {$j <= 2} {incr j} {
        set nodeTag [expr $i*3 + $j + 37]
        node $nodeTag [expr $i*$L1] [expr $j*$L2] [expr 3*$H]
    }
}

# 边界约束 - 底层全固定
for {set i 1} {$i <= 12} {incr i} {
    fix $i 1 1 1 1 1 1
}

# 几何变换（为不同方向的梁定义不同变换）
geomTransf Linear 1 1 0 0;  # 柱单元变换，局部y轴沿全局x轴
geomTransf Linear 2 0 1 0;  # X方向梁单元变换，局部z轴沿全局y轴
geomTransf Linear 3 0 0 1;  # Y方向梁单元变换，局部z轴沿全局z轴

# 截面几何参数（严格按文档H型钢尺寸计算，国际单位制：米）
# 梁截面 H200×600×20×20（翼缘宽×梁高×翼缘厚×腹板厚）
set bf_beam 0.200;      # 翼缘宽度 m
set h_beam 0.600;       # 梁高 m
set tf_beam 0.020;      # 翼缘厚度 m
set tw_beam 0.020;      # 腹板厚度 m
# 梁截面特性计算（精确计算）
set A_beam [expr 2*$bf_beam*$tf_beam + ($h_beam-2*$tf_beam)*$tw_beam];  # 截面面积 m²
set Iz_beam [expr ($bf_beam*pow($h_beam,3) - ($bf_beam-$tw_beam)*pow($h_beam-2*$tf_beam,3))/12];  # 强轴惯性矩 m⁴
set Iy_beam [expr (2*$tf_beam*pow($bf_beam,3) + ($h_beam-2*$tf_beam)*pow($tw_beam,3))/12];  # 弱轴惯性矩 m⁴
set J_beam [expr 2*$bf_beam*pow($tf_beam,3)/3 + ($h_beam-2*$tf_beam)*pow($tw_beam,3)/3];  # 扭转常数 m⁴

# 柱截面 H600×600×20×20（翼缘宽×柱高×翼缘厚×腹板厚）
set bf_col 0.600;       # 翼缘宽度 m
set h_col 0.600;        # 柱高 m
set tf_col 0.020;       # 翼缘厚度 m
set tw_col 0.020;       # 腹板厚度 m
# 柱截面特性计算（精确计算）
set A_col [expr 2*$bf_col*$tf_col + ($h_col-2*$tf_col)*$tw_col];  # 截面面积 m²
set Iz_col [expr ($bf_col*pow($h_col,3) - ($bf_col-$tw_col)*pow($h_col-2*$tf_col,3))/12];  # 强轴惯性矩 m⁴
set Iy_col [expr (2*$tf_col*pow($bf_col,3) + ($h_col-2*$tf_col)*pow($tw_col,3))/12];  # 弱轴惯性矩 m⁴
set J_col [expr 2*$bf_col*pow($tf_col,3)/3 + ($h_col-2*$tf_col)*pow($tw_col,3)/3];  # 扭转常数 m⁴

# 材料常数（国际单位制：N/m²）
set E_steel 206.0e9;    # 弹性模量 N/m²
set G_steel 79.0e9;     # 剪切模量 N/m²

puts "按文档精确计算的截面特性（国际单位制：m）："
puts "梁截面H200×600×20×20: A=[format %.6f $A_beam]m², Iz=[format %.3e $Iz_beam]m⁴, Iy=[format %.3e $Iy_beam]m⁴"
puts "柱截面H600×600×20×20: A=[format %.6f $A_col]m², Iz=[format %.3e $Iz_col]m⁴, Iy=[format %.3e $Iy_col]m⁴"

# 创建纤维截面定义
# 柱截面（H600×600×20×20）- 钢-混凝土组合截面
section fiberSec 1 -GJ [expr $G_steel*$J_col] {
    # 混凝土核心区
    set y1 [expr -$h_col/2 + $tf_col]
    set z1 [expr -$bf_col/2 + $tw_col]
    set y2 [expr $h_col/2 - $tf_col]
    set z2 [expr $bf_col/2 - $tw_col]
    patch rect 1 10 10 $y1 $z1 $y2 $z2
    
    # 钢材翼缘（上下）
    patch rect 2 10 2 [expr -$h_col/2] [expr -$bf_col/2] [expr -$h_col/2 + $tf_col] [expr $bf_col/2]
    patch rect 2 10 2 [expr $h_col/2 - $tf_col] [expr -$bf_col/2] [expr $h_col/2] [expr $bf_col/2]
    
    # 钢材腹板（左右）
    patch rect 2 2 10 [expr -$h_col/2 + $tf_col] [expr -$bf_col/2] [expr $h_col/2 - $tf_col] [expr -$bf_col/2 + $tw_col]
    patch rect 2 2 10 [expr -$h_col/2 + $tf_col] [expr $bf_col/2 - $tw_col] [expr $h_col/2 - $tf_col] [expr $bf_col/2]
}

# 梁截面（H200×600×20×20）- 纯钢截面
section fiberSec 2 -GJ [expr $G_steel*$J_beam] {
    # 钢材翼缘（上下）
    patch rect 2 10 2 [expr -$h_beam/2] [expr -$bf_beam/2] [expr -$h_beam/2 + $tf_beam] [expr $bf_beam/2]
    patch rect 2 10 2 [expr $h_beam/2 - $tf_beam] [expr -$bf_beam/2] [expr $h_beam/2] [expr $bf_beam/2]
    
    # 钢材腹板
    patch rect 2 2 20 [expr -$h_beam/2 + $tf_beam] [expr -$tw_beam/2] [expr $h_beam/2 - $tf_beam] [expr $tw_beam/2]
}

# 积分点数定义
set numIntgrPts 5;

# 柱单元定义
# 1-2层柱
for {set i 1} {$i <= 12} {incr i} {
    element forceBeamColumn [expr 100+$i] $i [expr $i+12] 1 Lobatto 1 $numIntgrPts
}

# 2-3层柱
for {set i 13} {$i <= 24} {incr i} {
    element forceBeamColumn [expr 100+$i] $i [expr $i+12] 1 Lobatto 1 $numIntgrPts
}

# 3-4层柱
for {set i 25} {$i <= 36} {incr i} {
    element forceBeamColumn [expr 100+$i] $i [expr $i+12] 1 Lobatto 1 $numIntgrPts
}

# 梁单元定义
# 二层梁
# X向梁（沿X方向，使用变换2）
for {set j 0} {$j <= 2} {incr j} {
    for {set i 0} {$i <= 2} {incr i} {
        set node1 [expr 13 + $i*3 + $j]
        set node2 [expr 13 + ($i+1)*3 + $j]
        set elemTag [expr 200 + $i*3 + $j + 1]
        element forceBeamColumn $elemTag $node1 $node2 2 Lobatto 2 $numIntgrPts
    }
}

# Y向梁（沿Y方向，使用变换3）  
for {set i 0} {$i <= 3} {incr i} {
    set node1 [expr 13 + $i*3]
    set node2 [expr 13 + $i*3 + 1]
    set elemTag [expr 210 + $i + 1]
    element forceBeamColumn $elemTag $node1 $node2 3 Lobatto 2 $numIntgrPts
    
    set node1 [expr 13 + $i*3 + 1]
    set node2 [expr 13 + $i*3 + 2]
    set elemTag [expr 215 + $i + 1]
    element forceBeamColumn $elemTag $node1 $node2 3 Lobatto 2 $numIntgrPts
}

# 三层梁（复制二层梁，节点编号+12）
# X向梁
for {set j 0} {$j <= 2} {incr j} {
    for {set i 0} {$i <= 2} {incr i} {
        set node1 [expr 25 + $i*3 + $j]
        set node2 [expr 25 + ($i+1)*3 + $j]
        set elemTag [expr 300 + $i*3 + $j + 1]
        element forceBeamColumn $elemTag $node1 $node2 2 Lobatto 2 $numIntgrPts
    }
}

# Y向梁
for {set i 0} {$i <= 3} {incr i} {
    set node1 [expr 25 + $i*3]
    set node2 [expr 25 + $i*3 + 1]
    set elemTag [expr 310 + $i + 1]
    element forceBeamColumn $elemTag $node1 $node2 3 Lobatto 2 $numIntgrPts
    
    set node1 [expr 25 + $i*3 + 1]
    set node2 [expr 25 + $i*3 + 2]
    set elemTag [expr 315 + $i + 1]
    element forceBeamColumn $elemTag $node1 $node2 3 Lobatto 2 $numIntgrPts
}

# 四层梁（复制二层梁，节点编号+24）
# X向梁
for {set j 0} {$j <= 2} {incr j} {
    for {set i 0} {$i <= 2} {incr i} {
        set node1 [expr 37 + $i*3 + $j]
        set node2 [expr 37 + ($i+1)*3 + $j]
        set elemTag [expr 400 + $i*3 + $j + 1]
        element forceBeamColumn $elemTag $node1 $node2 2 Lobatto 2 $numIntgrPts
    }
}

# Y向梁
for {set i 0} {$i <= 3} {incr i} {
    set node1 [expr 37 + $i*3]
    set node2 [expr 37 + $i*3 + 1]
    set elemTag [expr 410 + $i + 1]
    element forceBeamColumn $elemTag $node1 $node2 3 Lobatto 2 $numIntgrPts
    
    set node1 [expr 37 + $i*3 + 1]
    set node2 [expr 37 + $i*3 + 2]
    set elemTag [expr 415 + $i + 1]
    element forceBeamColumn $elemTag $node1 $node2 3 Lobatto 2 $numIntgrPts
}

# 质量定义（按文档楼板荷载计算，国际单位制：kg）
# 楼板组合荷载：6.5kN/m² = 6500N/m²
# 每个节点影响面积：4.0m × 3.5m / 4 = 3.5m²（平均分配给周围4个节点）
# 节点重量：6500N/m² × 3.5m² = 22750N
# 节点质量：22750N / 9.8m/s² = 2321.4kg
set nodeMass 2321.4
# 二层节点 (13-24)
for {set i 13} {$i <= 24} {incr i} {
    mass $i $nodeMass $nodeMass $nodeMass 0 0 0
}
# 三层节点 (25-36)  
for {set i 25} {$i <= 36} {incr i} {
    mass $i $nodeMass $nodeMass $nodeMass 0 0 0
}
# 四层节点 (37-48)
for {set i 37} {$i <= 48} {incr i} {
    mass $i $nodeMass $nodeMass $nodeMass 0 0 0
}

# 静力荷载（按文档组合荷载6.5kN/m²计算，国际单位制：牛顿）
# 每个节点荷载 = 组合荷载 × 影响面积
# 影响面积 = 4.0m × 3.5m / 4 = 3.5m²
# 节点荷载 = 6.5kN/m² × 3.5m² = 22.75kN = 22750N
set nodeLoad 22750.0

# 测试1：Z向荷载（重力方向）
pattern Plain 1 Linear {
    # 二层节点荷载 (13-24)
    for {set i 13} {$i <= 24} {incr i} {
        load $i 0 0 [expr -$nodeLoad] 0 0 0
    }
    # 三层节点荷载 (25-36)
    for {set i 25} {$i <= 36} {incr i} {
        load $i 0 0 [expr -$nodeLoad] 0 0 0
    }
    # 四层节点荷载 (37-48)
    for {set i 37} {$i <= 48} {incr i} {
        load $i 0 0 [expr -$nodeLoad] 0 0 0
    }
}

# 静力分析
constraints Transformation
numberer RCM
system BandGeneral
test NormDispIncr 1e-6 100 0 1
algorithm Newton
integrator LoadControl 0.1
analysis Static

puts "开始静力分析..."
analyze 10
loadConst -time 0.0
puts "静力分析完成"

# 记录关键节点位移（国际单位制：米）
puts "\n静力分析结果："
puts "节点17（二层②轴）: Z位移 = [format %.6f [nodeDisp 17 3]] m"
puts "节点29（三层②轴）: Z位移 = [format %.6f [nodeDisp 29 3]] m" 
puts "节点41（四层②轴）: Z位移 = [format %.6f [nodeDisp 41 3]] m"

# 检查所有楼层节点的位移
puts "\n各楼层位移检查（单位：米）："
puts "二层所有节点Z位移："
for {set i 13} {$i <= 24} {incr i} {
    puts "节点$i: Z=[format %.6f [nodeDisp $i 3]] m"
}
puts "\n三层所有节点Z位移："
for {set i 25} {$i <= 36} {incr i} {
    puts "节点$i: Z=[format %.6f [nodeDisp $i 3]] m"
}

# 地震动分析设置
wipeAnalysis
remove loadPattern 1

# 地震动时程设置（单位：m/s²）
set dt 0.001;  # 时间步长 0.01s

# 创建时间序列（使用elcentro地震动，参考算例缩放因子）
timeSeries Path 1 -dt $dt -filePath "earthquake.txt" -factor 150.0

# 创建地震动荷载模式（X方向地震动）
pattern UniformExcitation 2 1 -accel 1

# 瑞利阻尼设置（按模态阻尼比5%）
set eigenValues [eigen 2]
set lambda1 [lindex $eigenValues 0]
set lambda2 [lindex $eigenValues 1]
set omega1 [expr sqrt($lambda1)]
set omega2 [expr sqrt($lambda2)]
set xi 0.05;  # 阻尼比5%
set a0 [expr $xi*2.0*$omega1*$omega2/($omega1+$omega2)]
set a1 [expr $xi*2.0/($omega1+$omega2)]
rayleigh $a0 0.0 0.0 $a1

# 记录水平位移时程（X方向）
recorder Node -file "displacement_x.out" -time -node 17 29 41 -dof 1 disp

# 记录水平加速度时程（X方向）
recorder Node -file "acceleration_x.out" -time -node 17 29 41 -dof 1 accel

# 动力分析设置
constraints Transformation
numberer RCM
system BandGeneral
test NormDispIncr 1e-6 25 0
algorithm ExpressNewton 2 1 1.0
integrator Newmark 0.55 0.275625
analysis Transient

# 执行地震分析（按参考算例方式）
set analysisTime 40.0
set numSteps [expr int($analysisTime/$dt)]
puts "开始地震动时程分析，共$numSteps步..."
set dynamicStartTime [clock milliseconds]
analyze $numSteps $dt
set dynamicEndTime [clock milliseconds]
set dynamicTime [expr $dynamicEndTime - $dynamicStartTime]

puts "地震动时程分析完成"
puts "动力分析耗时: ${dynamicTime} 毫秒 ([format %.3f [expr $dynamicTime/1000.0]] 秒)"

# 输出最终位移结果（单位：米）
puts "\n地震动分析最终位移结果（单位：米）："
puts "节点17（二层②轴）: X位移 = [format %.6f [nodeDisp 17 1]] m"
puts "节点29（三层②轴）: X位移 = [format %.6f [nodeDisp 29 1]] m"  
puts "节点41（四层②轴）: X位移 = [format %.6f [nodeDisp 41 1]] m"

# 模态分析
set eigenValues [eigen 5]
puts "\n模态分析结果："
for {set i 0} {$i < [llength $eigenValues]} {incr i} {
    set eigenVal [lindex $eigenValues $i]
    if {$eigenVal > 0} {
        set freq [expr sqrt($eigenVal)/(2*3.14159)]
        set period [expr 1.0/$freq]
        puts [format "第%d阶: 频率 = %.3f Hz, 周期 = %.3f 秒" [expr $i+1] $freq $period]
    }
}

puts "\n模型创建完成！（国际单位制：m-N-kg）"